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I. INTRODUCTION 


The work performed by the author while at the Langley Research 
Center (LRC) during the period June 1, 1977 through May 31, 1978 on 
an IPA assigment agreement provided the necessary continuity between 
the work carried out under the grant NSG 47-004-114 and the present 
grant NSG 1546. Several ideas evolved or were stimulated as a result 
of author's collaboration with the technical monitors Robert Hayduk and 
Robert Thomson of LRC. 

A multitude of problems common to the development of a successful 
nonlinear analysis code prevented ACTION. from being able to solve even 
the simplest of the nonlinear problems during the early stages. Most 
of these problems having been overcome ACTION'S performance looked cau- 
tiously optimistic. Hence, the need to validate the ACTION computer 
code on an aircraft-like structure followed by its finalization for pub- 
lic release were recognized. 
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II. TECHNICAL PROGRESS 


The technical efforts during the past year were concentrated in 
several different areas. These are discussed below. Each group of 
paragraphs provides a background, a brief statement of the problem, 
approach to solution, the progress made and is some cases recommenda- 
tions for future work. 

1. Structuring of the ACTION Code to Handle Relatively Large Scale 
Problems 

In order to make ACTION suitable to handle large scale problems 
involving finite element models with stringer, frame and membrane ele- 
ments, it was necessary to undertake the following tasks. 

a. Implementation of Analytic Gradients for the Membrane Element 
The calculation of gradients of the strain energy of deformation of 
the membrane element was hitherto being performed by the use of central 
differences in ACTION. This is nearly twice as more expensive as when 
the gradients are calculated analytically. This is especially crucial . 
when a model contains a large number of membrane elements. It was thus 
imperative that gradients of strain energy of deformation of the mem- 
brane elements be calculated analytically in ACTION. Accordingly, on 
lines similar to those used for the frame element explicit expressions 
were developed for the calculation of the gradients of the strain energy 
both in the elastic and the inelastic range. In the inelastic range the 
computations were simplified by the use of the total deformation theory 
[1]. If incremental flow theory of plasticity is to be used then the 
definition of the potential function must be altered to allow the use of 

I 

certain minimum principles in plasticity '[2]. This novel formulation 
using the direct minimization approach appears like a topic suitable for 
Master's thesis and hence Mr. Lin T. Duong who has been working as a 
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graduate research assistant on the project indicated his interest in 
pursuing this line of research towards his Master of Science degree. 
Whether such a formulation will be computationally more efficient or 
will result in improved response prediction remains to be seen. 

b. Quasi-Newton Minimization Algorithms that Exploit Sparsity 
Of the two minimization algorithms available in ACTION, Powell's 
conjugate algorithm can handle problems involving thousands of variables 
but its convergence rate is at best linear. Hence, it is unlikely that 
it can be competitive with second order algorithms like the Fletcher's 
algorithms . (BFGS) [3] whose convergence, rate is at least superlinear. 
Fletcher's method is especially suited for solution of transient pro- 
blems in steps since it has a reasonably good approximation to the 
variables and the inverse Hessian of the potential surface. The number 
of minimizations required for convergence to a solution after the first 
time step is thus a very small fraction of the total number of degrees 
of freedom. The drawback of Fletcher's algorithm is however it's stor- 
age requirements - N*(N+l)/2 for a N degree of freedom problem. This 
is because Fletcher's algorithm approximates the Hessian inverse rather 
than the Hessian itself which is usually banded and very sparse. Thus, 
if an extremely large scale nonlinear transient problem is to be solved 
using energy minimization technique it is imperative to use a variable 
metric algorithm which updates the Hessian rather than its inverse and - 
thereby exploits and maintains its sparsity in its march, to the minimum. 
It appeared at one point that for solving problems with as few as 325 de- 
grees of freedom using Fletcher's method the available core on the CDC 
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system may be insufficient.* Accordingly, work was initiated to im- 
plement into ACTION a minimization algorithm that exploits sparsity. 

Such’ an algorithm based on the work of Sphubert [4] has been implemented 
into one of VPI & SU’s version of ACTION. The algorithm, however, 
remains to be validated. A survey of the state-of-the-art of minimiza- 
tion techniques for extremely large scale problems, which was conducted 
in this connection, may be found in reference [5]. 

2. Drop Test of the Navajo Fuselage Section : 

VPI & SU was entrusted with the task of development of a model for 
nonlinear crash dynamic analysis of the Navajo fuselage section using 
the ACTION code. The response characteristics of such a model was to be 
compared with comparable models developed using DYCAST [6] and KRASH 
[7]. Such a comparison was subsequently reported in reference [8]. 

An ACTION model of the Navajo fuselage section was developed using 
the corresponding DYCAST model (an initial cruder model different from 
the one reported in [8]) as its basis. This model shown in Fig. 1 has 
105 nodes, 209 members (36 stringers, 77 frame elements and 96 mem- 
branes) , 71 lumped masses and a total of about 336 degrees of freedom. 
This is the largest model ever attempted using ACTION. 

Although the model could be run using the Powell's conjugate 
gradient minimization • algorithm (MIN5) which requires very little stor- 
age, it was believed, as explained before, that in going from one step 

to the next this algorithm does not have the same advantages of the 

*An ACTION version which could handle problems with up to about 475 
degrees of freedom using the Fletcher's method was put together at VPI 
& SU. This was possible because of the almost unlimited storage avail- 
able on the IBM system. Although some preliminary runs were made using 
this version the task of simulation of large scale problems like the 
Navajo drop test had to be relegated to the CDC computer at LRC because 
of extremely tight funding situation at VPI & SU'. 



Drop test of eight passenger plane, fusetage section 

Action mode! 



second order algorithm of Fletcher (MIN7). The latter has a good ap- 
proximation to the inverse Hessian in addition to a good initial guess 
for the variables at hand at the end of each time step. 

To run this model using the more efficient Fletcher's minimization 
algorithm in ACTION it was necessary to overlay the ACTION code on 
NASA's CDC system using the segloader. The assistance of Ms, Barbara 
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Durling of NASA Langley in accomplishing this task is greatly appreciated. 

Of simulation interest was the occupant chest motion and its ver- 
tical acceleration at the pelvis location. The occupant was modeled by 
a single lumped mass while the seat was modeled by a set of four non- 
linear stringer elements whose stress-strain behavior was based on 
previous, independent static crash tests on similar seats. This infor- 
mation was provided by LRC. Although, the ground plane capability 
(terrain model) in ACTION could have been exercised for this problem, in 
the interest of simplicity, the aircraft substructure was assumed to be 
in contact with the. ground plane at nodes (86) and (91) while a velocity 
of 330 inches/sec. was imparted to the entire model. Thus one would 
expect correlation only in' the initial phases of the response. 

Figures 2 through 5 provide the correlation between the analysis 
and test results. For additional model and simulation details it is 
suggested that reference [8] be consulted. This reference also provides 
a comparison of the performance of the energy minimization technique 
vis-a-vis the so-called hybrid technique (KRASH) and another technique 
which utilizes the pseudo force technique (DYCAST) . 
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Fig. 3 Comparison of occupant Pelvis vertical 
Acceleration , Test versus Analysis 


8 





Fig. 4- Comparison of fuselage floor vertical 
acceleration , Test versus Analys/s 
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Fig. 5 Comparison- of fuselage floor vertical 
acceleration } Test versus Analysis. 
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3. Implementation of the Frame Element .with Cross-Sectional Flexi- 


bility into DYCAST 

VPI & SU was entrusted with the task of aiding Greenman Engineers 
in implementing the frame element with cross-sectional flexibility into 
their program DYCAST. The frame element with cross-sectional flexibility 
was developed by VPI &. SU in the ACTION code environment [9] which uses 
the direct minimization approach. Implementation of this element re- 
quires the development of a tangent stiffness matrix for such an ele- 
ment. This task was given a very low priority dictated by the conclu- 
sions of the study ■ reported in reference [ 9 ] and other experiences in 
simulating the crash response of the Navajo section. 

Report [9] concludes that the frame elements in both ACTION and 
DYCAST which do not currently allow for cross-sectional deformations are 
adequate in predicting gross response parameters even though the struc- 
tures in question may undergo severe localized cross-sectional deforma- 
tions. By gross parameters is meant total energies, load-deflection 
response, etc. The report goes on to conclude that purely from the 
point of view of crashworthiness where the trauma measures are deter- 
mined more by gross parameters it seems highly unlikely that they will 
be significantly influenced by highly localized cross-sectional deforma- 
tions. The inclusion of such effects can only make the already expen- 
sive nonlinear analysis only more so without any. significant pay-off by 
way of improvement in response. Considering the cost comparisons of 
KRASH, ACTION & DYCAST in reference [8] vis-a-vis the quality of their 
respective response predictions this would warrant a careful assessment 
of the necessity of including such a costly element into either ACTION 
or DYCAST. In any event if such am implementation is to be carried out 
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into DYCAST at some future date, VPI & SU’s cooperation, in this matter 
can always be counted upon. The formulation details outlined in ref- 
erence [9] will provide the required basis for this purpose. 

4.. Documentation and Wrap-up of the ACTION Code 

a. ACTION Theory Document 

The old ACTION theory 'document was thoroughly revised to reflect 
the current formulation basis in ACTION (especially the part pertaining 
to the calculation of analytic gradients and the new minimization algo- 
rithms) . A copy of the completed- rough draft of this document will be 
mailed to LRC within the next few days. Upon receipt. of comments from 
LRC the document will be revised accordingly and the final version of 
this document will be completed under NASA's Master Agreement arrange- 
ment with VPI & SU (NAS1-15080 Task #10). 

b. , ACTION User's Guide 

The existing ACTION user's manual was revised to be consistent with 
the final ACTION version to be released through COSMIC. Additional 
features included analytic gradients, selection of Newmark-Beta para- 
meters, automatic calculation of lumped masses, a new solid rectangular 
frame cross-section, etc- Several outdated, unvalidated ACTION features 
have been removed from the present version as called for in the grant 
proposal. The revised ACTION version, however, remains to be thoroughly 
tested. This testing will be performed again under the Master Agreement 
arrangement contract NAS1-15080 Task #10. A copy of the completed rough 
draft of this document will be mailed to LRC within the next few days 
for their comments to be implemented in the final version. 
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c. Publication, of VPI&SU Reports and Papers in Referred Journals . 

Progress in this area was extremely satisfactory. A report [9] en- 
titled "An Investigation into the Effects of Cross-Sectional Flexibility 
on Response" VPI-E-78-30 was completed and four copies of the same were 
mailed to LRC (see Appendix A) . This study summarizes the findings 
of the importance of cross-sectional flexibility on gross response and 
attempts to identify reasons for the discrepancy between analytical 
and experimental predictions. A technical note entitled "On the Impor- 
tance of Cross-Sectional Flexibility on Gross Response" has been ac- 
cepted for publication by the Journal of Computers and Structures (see 
Appendix B) . 

A synoptic summarizing the ACTION formulation and validation en- 
titled "Nonlinear Transient Analysis via Energy Minimization” has been 
accepted for publication by the Journal of American Institute of Aero- 
nautics and Astronautics. The full length paper by the same title has 
been accepted as a backup NTIS document number N79-22819 for the synop- 
tic. The same full length L was previously released as a VPI & SU report 
number VPI-E-79-10 entitled "Nonlinear Transient Analysis of Aircraft- 
like Structures - Theory and Validation" (see Appendix C) . The nec- 
essary number of copies were mailed to LRC. 

A paper entitled "Energy Minimization versus Pseudo Force Technique 
for Nonlinear Structural Analysis" has been accepted for publication 
by the Journal of Computers and Structures (see Appendix D) . This paper 
compares the efficiency of the minimization technique (ACTION) vis-a-vis 
the pseudo force technique (DYCAST) for the solution of nonlinear pro- 
blems. It also explores the feasibility of extending the minimization 
technique for the solution of extremely large scale nonlinear problems 
competitively. 



5. Simplified Modeling Techniques 

This is the only task which received the least attention because of 
lack of time and also because of lack of suitable personnel who could be 
assigned to tackle this problem. Some nominal progress was, however, 
made. A literature survey on the topic revealed that ideas similar to 
those proposed in the proposal have been explored in much greater detail 
by Kawai and his associates [10]. Their work suggests that ideas like 
the ones proposed in the proposal should be worthwhile to pursue es- 
pecially if codes like DYCAST or ACTION can ever hope to be cost effec- 
tive with codes like KRASH [8]. 
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III. OPERATIONAL SET-UP 


The work reported in this report was carried out by the principal 
investigator, Dr. M. P. Kamat, with the aid of Mr. Linh T; Duong, a 
graduate research assistant in the department of Engineering Science and 
Mechanics of VPI & SU- Because of lack of students with the appropriate 
background and the time required to familiarize them with the ACTION 
code no attempt was made to hire other graduate research assistants. 
Instead the principal investigator undertook a larger share of the 
responsibility inorder that he may meet the grant obligations more 
fully. 
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IV. CONCLUSIONS 


The work performed under this grant achieved several objectives. 

It successfully demonstrated ACTION’S capability to model and analyze 
aircraft for nonlinear crash response at a cost which is comparable 
with that of other similar analytical tools like DYCAST. It was estab- 
lished that computer programs like ACTION & DYCAST in their present 
status are quite adequate in predicting gross response parameters in- 
spite of the fact that structures in question may undergo severe local- 
ized cross-sectional deformations. The theoretical and user’s documen- 
tation for the ACTION code were nearly finalized for public release. In 
addition to a few technical reports several papers were published in 
referred Journals. State-of-the-art survey conducted indicates minimi- 
zation techniques to be quite suitable for the solution of large scale 
nonlinear px*oblems. A similar survey also indicates that simplified 
modeling techniques are not only viable but should be pursued if pro- 
grams like ACTION or DYCAST are to be cost effective with programs like 
KRASH. 
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Abstract 


A simplified' frame element (3D-beam) formulation which accounts 
for the beam cross-sectional flexibility is presented. The effectiveness 
of such a formulation as opposed to modeling the frame element using 
membranes and plate bending elements is investigated on the post-buckling 
response of a beam column with a thin-walled channel cross-section. 

The fidelity of the simulation with and without such effects is determined 
on the basis of a comparison with available experimental results and 
other -independent simulations . The study concludes that in view of the 
excellent correlation between experimental and mathematical predictions 
using simulators like ACTION and O-PLANE-MG which do not account for 
cross-sectional flexibility, a beam element with cross-sectional 
flexibility will have little. pay-off , especially so for the response 
parameters of importance in crashworthiness.. 
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1 . Introduction 


Static and dynamic load tests on tubular and angular frame structures 
conducted by the Dynamic Loads Branch of the NASA Langley Research Center 
reveal that significant cross-sectional warping and distortions occur 
near the joints [1]. It is believed that the lack of correlation between 
experimental and theoretical predictions may be attributed, at least in 
part, to cross-sectional flexibility during deformation hitherto not 
.accounted for by most mathematical simulators of such nonlinear phenomena. 
The present investigation was thus initiated with the objective of 
implementing such a capability into existing programs for nonlinear 
analysis of structures. 

In the context of the cross-sectional flexibility the word "beam" is 
perhaps a misnomer since it is implicit in its usage that its response is 
predominantly in the direction of its longest dimension which is supposedly 
an order of magnitude larger than its other cross-sectional dimensions. 
Hence when cross-sectional flexibility effects are significant it would 
be more appropriate to refer to the structure as either a three dimensional 
solid or in the case of thin-walled open or closed sections as an 
assemblage of flat or curved plates. In other words, cross-sectional 
flexibility in the case of thin-walled structures would be most naturally 
accounted for by idealizing the member or structure with membranes and 
plate bending finite elements. However, the degrees of freedom of the 
resulting model would in most cases, be prohibitive for a general nonlinear 
analysis. Hence, some other cheaper formulation which perhaps sacrifices 
some of the rigor of the formulation would be desirable. 



A thin -walled open section frame (3D-beam) element formulation which 
abandons the assumptions of the elementary beam theory and accounts for 
the flexibility of its cross-section in an approximate manner was undertaken. 
The objective was to develop a model which is simple and in terms of its 
size well below that which would result from building up the frame element 
using membrane and plate bending elements. Consequently, certain simplifying 
assumptions were incorporated in its development. The development was 
initially carried out in the context of the ACTION simulator. This 
simulator uses the minimization of the total potential energy for 
nonlinear analysis [2] . , 
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2. Formulation 


a. Assumptions : Restricted to beams with thin-walled cross-sections, 

the formulation focusses on distortions in the plane of the cross-sections 
only since warping, at least of the unrestrained type is usually accounted 
-for in most conventional frame element development. This. is true of the 
ACTION simulator as well. 

The development consists in abandoning the assumption of the 
elementary beam theory about cross-sections remaining rigid during 
deformation. Co-ordinates of a select few reference points on the cross- 
section are treated as additional degrees of freedom with a prescribed 
co-ordinate variation in between reference points. In essence then this 
is tantamount to treating the co-ordinates of the integration or 
quadrature points used for the evaluation of the strain energy of 
deformation of the frame element as additional degrees of freedom of the 
model. Since the problem now becomes extremely nonlinear it is reasonable 
to adopt a linearization procedure of the following type. At the beginning 
of each load increment elementary beam theory is used to predict overall 
deformations of the frame element assuming a rigid cross-section. At the 
end of the increment the load is maintained constant and the additional 
deformations of the elements of the cross-section are determined by 
treating the elements as prestressed plate elements with a known initial 
eccentricity. The shape of the cross-section is determined by minimizing 
the corresponding potential energy of additional deformations which consists 
of the strain energy of additional deformations and the potential of the 
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prestress. Using this new cross-section the process is then repeated for 
the next load increment. Iteration at constant, load could be employed to 
obtain improved results. 

b. Formulation Details : Consider the IE cross-section shown in 

Figure 1. This cross-section has five flange pieces and two web pieces. 

This element as implemented in the ‘ACTION Simulator has the option of 
specifying zero thicknesses for all but the two web pieces. This enables 
the generation of cross-sections of many well-known shapes. Each of the 
seven pieces is treated like a flat plate element extending between the 
two end nodes of the beam or frame element. A total of seven reference 
points^!) through(7)as shown in Figure 1 are chosen on the cross-section 
at each end of the beam element and half-way between the two nodes. 
Incidently, these locations of reference points are coincident with the 
Lobatto integration points (at the ends of the intervals of integrations) 
used for each of the seven plate elements. The stress-strain history at 
all of these points is thus available without additional calculation 
effort or storage. The cross-section is allowed to deform relative to 
the fixed point © in the cross-section. At each of the reference 
points ©, ©> ©> © and © only one independent displacement 
transverse to the respective plate elements CD , CD , 0 , ED and Q is 
allowed and at each of the reference points © and © only one independent 
displacement transverse to the respective plate elements JH and HI is 
allowed. The other position co-ordinate of each of the reference points 
is determined on the basis of the additional deformations of the plate 
elements being, inextensional . 
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For the IE section of Figure I the potential energy .of aciditionai 
deformations, it, is thus given by [3] 

ir = iTi + ir 2 (1-a] 

where 
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The deformation shapes w^ and' v^ are chosen as in the finite element 

formulation namely as polynomials with the values of the independent 

co-ordinates at the six reference points as coefficients of the 

interpolants. Thus for a typical element shown in Figure 2-a 
6 
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6 . ■ 

■and v. = 7 v..N.f£,p) i = 3 and 5 

i L-i ij } 

3=1 


(2-a) 


(2-b) 


- 5 - 



where ’s are the interpolating polynomials given below 

n : , = h (1-y) C? 2 -?) • 
n 2 = % (1+y)(? 2 -?) 
n 3 = % Oy) Cl-? 2 ) 

N 4 = % (i+y) ( l-? 2 ) 

N s = % Cl-Y) (? 2 +?) 

Ng = % (1+Y) (? 2+ ?). 


(2-c) 
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i i 

i =- 1, 2, 4, 6 and 7 or i = 3 and 5. In evaluating in Equation (l _c ) 
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with v = 0.5. E^^ is the average tangent modulus for the ith plate 

element. This is taken to be the average of the tangent moduli >at all 
the nine Lobatto integration points for the element. This approximation 
enables an explicit evaluation of the integrals in the expression for 
Thus typically for the element shown in Figure 2-a 
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For i=3 or 5, w in Equation (4) is replaced by v. No such explicit 

expression is however possible for because of the occurrence of 

stress resultant terms N° . , N° . and N° . within the integrands. 

xxi xyi xzi 6 

it^ is therefore evaluated numerically using Lobatto quadratures. For 
three points in each direction Lobatto quadratures reduces to the 
well-known Simpson’s rule. 

c. Simplification of the model : The model described previously 

may be simplified by reducing the number of reference points per plate 

element from six to four. For a typical element shown in Figure 2~b 

equations (2) are then replaced by 
4 
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the interpolating polynomials given below 
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N 4 * h (1+0 (l+n) 

d. Limitations and Deficiencies of the Model : Firstly, like in 

most other simulators, the frame element material model 'in ACTION is 
strictly uniaxial, i.e. interaction between shear and normal stresses is 
ignored. Shear stresses are thus evaluated by using the shear flow 
theory of thin-walled structures assuming linear elastic behavior. 
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Secondly, with changes in cross-section the shift of the position 
of the shear center should be accounted for. The transverse forces 
and V 7 which are assumed to pass through the shear center .of the original 
rigid cross-section will then give rise to additional twisting moments. 
Because the calculation of the shear center for a thin-wal-led open 
cross-section of arbitrary shape is highly complicated in the interest 
of simplicity and efficiency this is avoided. 

Thirdly, f.or the assumed deformation patterns of the flanges and 
webs it is not possible to maintain the right angle bends between the 
flanges and webs at every point along the length. The deformation patterns 
imply that the plate elements are effectively hinged at some points along . 
intersections. To obtain, any sort of slope continuity it will be 
necessary to assume higher order interpolation functions in the n direction. 
This defeats the original objective of this model. Hence again in the 
interest of simplicity such a deficiency is tolerated. 
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3. Results and Conclusions 


a. Discussion of Results : 

For validating the proposed model the tests conducted by Mclvor et al 
on thin-walled open section beam columns and reported in reference [4] 
seemed appropriate. One such test involved the post-buckling response 
of a beam-column with a thin-walled channel cross-section as shown in 
Figure 3. To prevent failure by direct compression the column was tested 
at an inclination of 5° from the vertical with both ends of the column 
being clamped. A column under these conditions is extremely imperfection 
sensitive. Hence Mclvor suggests assuming an additional 1° offset for 
the mathematical model to simulate the inherent imperfections of the 
actual column tested. The post-buckling response of this column as 
observed in the experiment and as reported in [4] shows extremely severe 
cross-sectional deformations and is presumably a good problem for validation 
of the new beam element proposed herein. 

Before validating the proposed new beam element it was necessary to 
assess the importance of the magnitude of the degradation of response 
without allowing for such effects. This effort was initially undertaken. 

The results of this study are shown in Figures 3 and 4. Surprisingly 
enough the responses predicted by both ACTION and O-PLANE-MG [5] correlate 
extremely well with the experimental prediction. It would appear then 
that gross responses like load-deflection are not likely to be affected by 
localized cross-sectional deformations, irrespective. of their severity,, 
especially when the beam material models in both ACTION and O-PLANE-MG 
permit simulation of plastic hinges. This leads very naturally to the 
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investigation of the original problem, namely the NASA angular frame 
wherein large discrepancies between the mathematical and experimental 
predictions were responsible for initiating the study of the effects of 
beam cross-sectional flexibility. 

Before we get into a discussion of this study however, it may be 
appropriate to discuss the results of the attempts of using the new beam 
model in predicting the post-buckling response of the' beam-column of 
Figure 3. These attempts were not very successful. The simplified 
beam model with j=4 (see Eqs. 2-d § 2-e) failed to produce any significant 
cross-sectional deformations and had to be abandoned in place of the 
more refined model with j=6 (see Eqs. 2-a § 2-b) . The new beam elements 
were employed only between nodes ©-©>©-©> ©-©> ®-0 
and ©-©• This model although instrumental in predicting rather 
severe distortions of the cross-section troubles in converging to an 
accurate solution for the prescribed load steps, presumably because of 
an increased degree of nonlinearity, prevented completion of this study. 
Drastic reductions in load steps may overcome this problem but only so 
at a very high computational cost. Limited results of this study are 
tabulated in Table 1. To say the least, /it is disturbing to see that 
initially the column has a tendency to stiffen leading to a higher 
axial load carrying capability. This does not seem impossible however, 
because it is quite likely that the plastic section modulus of the 
slightly deformed beam cross-section may be higher than that of the 
undeformed cross-section. With increasing cross-sectional deformations 
however, the trend is likely to be reversed and this does appear to be 
characteristic of the results of Table 1. The cross-sectional shapes' 
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at different locations' along the length of the beam are illustrated 
qualitatively in Figure 5. While in some places these shapes appear to 
be intuitively reasonable in other places they do not. This may be the 
result of the simplicity of the model along with its underlying assumptions. 
Undoubtedly these results warrant further improvement and an investigation 
into the accuracy of the results generated using the new beam element. 

Next, as regards NASA's angular frame the lack of correlation 
between mathematical simulation and experimental results, if any, at least 
in the linear range must be attributed to either (i) lack of inappropriate 
boundary conditions that is to say boundary conditions of the mathematical 
model different from those in the experiment; (ii) assumptions regarding 
rigidity of the bulkheads' too conservative; (iii) insufficiently refined 
model or inappropriate discretization; (iv) lack of inclusion of some 
important response features like shear deformation. 

The results of a study entailing the sensitivity of response (due 
to' a 10 lb total load) to variations in modeling parameters are outlined 
in Table 2. The maximum deflection of the model of class 4 in Table 2 
agrees closely with that of the experimental model. Shear deformation 
effects which are likely to be important for this frame were not included 
in this analysis, however. With the inclusion of such effects the mathematical 
model will undoubtedly be even more flexible than the experimental model. 

It would appear that the correlation between the mathematical and 
experimental models can be substantially improved by using the class 4 
model of Table 2 or refinements thereof and allowing for shear flexibility. 

The frame in question is very much like a vierendel truss wherein the 
sheaf forces give rise to secondary bending. 
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The results of Table 2 suggest that the discrepancy between experi- 
mental and mathematical simulations as in Figure 6 are most likely. the result 
of inappropriate modeling of the angular frame and not the result of 

some highly localized cross-sectional deformations as was anticipated, 
b. Conclusions : 

It is clear from this study that computer programs like ACTION and 

0-PLANE-MG in their present status are quite adequate in predicting 

gross response parameters in spite of the fact that the structures in 

question may undergo severe localized cross-sectional deformations. 

By gross parameters we mean total energies, load-deflection responses etc. 

which are the quantities directly solved for in the displacement formulation. 

They are not derived quantities like strains, stresses etc. which are 

derived through spatial differentiation. If derived quantities are of 

importance we are talking of pointwise correlation rather than a global 

or gross correlation. In such an event the displacement formulations 

of programs like ACTION or 0-PLANE-MG even with the inclusion of the 

effects of the cross-sectional flexibility (in the manner outlined in this 

report) would be inadequate. One will have to then resort to some sort 

of hybrid or mixed formulations for better response predictions. 

* 

Purely from the point of view of crashworthiness where the trauma 
measures are determined more by gross parameters it seems highly unlikely 
that they will be significantly inf luenced- by highly localized cross- 
sectional deformations. The inclusion of such effects can only make the 
already expensive nonlinear analysis only more so without any 
significant pay-off by way of improvement ' in response. 
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TABLE 1. 


EFFECTS OF CROSS-SECTIONAL FLEXIBILITY (INITIAL STAGES): 


VERTICAL 

DEFL. 

IN 

LOAD P 

(KIPS) 


(p ) WITHOUT 
1 CSF 

/p x WITH 
P 2 CSF 

AP - P r - P 2 

0,03 

10.228 

10,217 

0.011 

0.08 

8.6639 

8,7798 

• -0.1159 

0,130 

6.0146 

6,2641 

-0.2495 

0,150 

5.5161 . 

5.9241 ' 

-0.408 

0.200' 

4,6866 

4,7734 

-0.0868 

0,225 

4.3972' 

4.6097 ■ 

-0.2125’ 



TABLE 2. 


SENSITIVITY OF RESPONSE TO MODELLING OF NASA'S ANGULAR FRAME 


MAX DEFL. x 10 3 in. 

TYPE OF MODEL FOR P-10 lb, 


ACTION 


1. a, THREE D.O.F. PER NODE 7.5451 

b. BLKHDS. NOT RIGID 

c, END FULLY FIXED 


2, a. SIX D.O.F. PER NODE 8.3433 

b. BLKHDS. NOT RIGID 

c. END FULLY FIXED 


3. A. SIX D.O.F. PER NODE 5.6992 

b, BLKHDS. RIGID 

c. END FULLY FIXED 




4, a. SIX D.O.F,. PER NODE 16,5658 

b, BLKHDS. NOT RIGID 

c. END REPLACED BY A FLEXIBLE 
' BLKHD . WITH PROPER SUPPORT 


On the Importance of Cross-sectional 
Flexibility on Gross Response* 

M. P. Kamat 

Department of Engineering Science & Mechanics 
Virginia Polytechnic Institute and State University 
Blacksburg, Virginia 24061 

ABSTRACT 

Most nonlinear analyzers like ACTION and PLANS do not account 
for cross-sectional flexibility of thin-walled frame members. They 
do however, permit -simulation of plastic hinges. Judging by the 
excellent correlation between theory and experiment, such analyzers 
appear quite adequate in predicting reasonably accurately gross re- 
sponse parameters of thin-walled frame structures even though they 
may undergo severe cross-sectional deformations. 

INTRODUCTION . 

The crush response of structural components of an automobile or 
aircraft which are built-up from thin-walled open section frame ele- 
ments entails severe deformations of the cross-section. Static and 
dynamic load tests on tubular and angular frame structures conduct- 
ed by the Dynamic Loads Branch of the NASA Langley Research Center 
revealed that significant cross-sectional warping and distortions - 
occur near the joints [1 ] . In fact, it was conjectured that the 


* This work was supported by NASA Langley under grant NGR 47-004-114 
under the congnizance of the technical monitor, Robert J. Hayduk. 



large discrepancies between experimental and theoretical predictions 
may, at least in part, be attributable to such phenomena which are 
normally not accounted for in a conventional nonlinear finite ele- 
ment analysis. 

Recently, attempts have been made to model the phenomenon of 
cross-sectional deformations analytically [2], [3]. The purpose of 
this paper however, is not to evaluate these attempts but rather to 
assess the degradation of the fidelity of mathematical simulation of 
the response using nonlinear analyzers like ACTION [4] and PLANS [5] 
which currently do not permit such, model ing. 

RESULTS AND CONCLUSIONS 

Anderson, Mclvor and Kimball have conducted several tests on 
thin-walled open section columns and these have been reported in ref- 
erence [6]. These tests seem appropriate for the present investiga- 
tion. One such test involved the post-buck! ing elastic-plastic re- 
sponse of a beam-column with a thin-walled channel cross-section as 
shown in Figure 1. To prevent failure by direct compression, the 
column was tested at an inclination of 5° from the vertical with both 
ends of the column being clamped. A column under these conditions 
is extremely imperfection sensitive. Hence, Anderson et al [6] 
suggest assuming an additional 1° offset for the mathematical model 
to simulate the inherent imperfections of the actual column tested. 
The post-buckling response of this column as observed in the experi- 
ment and as reported in reference [6] involves extremely gross cross- 
sectional deformations and is presumably a good problem for the pur- 
poses of this paper. 
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The column of Figure 1 exhibits a bifurcation phenomenon with 
a highly unstable branch. For reasons very well known, most nonlin- 
ear finite element analyzers like ACTION and PLANS have to resort to 
a displacement rather than a load incrementation for predicting re- 
sponse along an unstable branch. The mathematical predictions of 
ACTION and PLANS have been illustrated in Figure 2 along with the 
experimental and that predicted by UMVCS-1 [7], 

Surprisingly enough, the responses predicted by ACTION and 
PLANS correlate extremely well with the experimental prediction. 

This excellent correlation between theory and experiment may be 
the result of several compensating assumptions, especially when one 
considers the fact that both ACTION and PLANS permit only large 
rotations with only small strains. The slight differences between 
the responses of ACTION and PLANS may be due to one or all of the 
following causes: (i) ACTION uses an energy minimization approach 

whereas PLANS uses a one-step Newton-Raphson psuedo force technique; 
(ii) the deformation models of ACTION and PLANS differ in the in- 
elastic range; (iii) the number of quadrature points for integration 
of energy densities and stresses are not identical. Itmust be re- 
marked, however, that with an elastic-perfectly plastic stress strain 
curve, both ACTION and PLANS do indeed possess the capabil ity of simu- 
lating a plastic hinge. UMVCS-1 is a program based on -the assump- 
tion that plastic deformations occur at idealized hinges located at 
nodal points. A generalized hinge theory which accounts for biaxial 
bending, .torsion and axial loading is employed. The user has to 
test the different types of joints and cross-sections experimentally 
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to gather data required for the constitutive laws employed in 

the analytical formulation. In this sense, it is more of a hybrid 

program. 

It is interesting to note that while the UMVCS-1 model is con- 
sistently stiffer than the actual that of PLANS is consistently more 
flexible than the actual and that of ACTION is more like an "average" 
model. If one were to account for cross-sectional flexibility in 
ACTION or PLANS, it would appear that this would lead to a deteriora- 
tion of the excellent correlation. 

It appears that the excellent agreement between theory and ex- 
periment in Figure 2, in spite of the absence of cross-sectional flex- 
ibility, is attributable to the capability of both ACTION and PLANS to 
permit simulation of plastic hinges. Gross responses like strain 
energy, load versus deflection are not likely to be affected by highly 
localized cross-sectional deformations, irrespective of their severity. 
It is the view of this author that nonlinear analyzers like ACTION 
and PLANS in their present status seem quite adequate in predicting 
reasonably accurately gross parameters like total strain energies, 
load-deflection responses, etc. of structures that undergo severe 
localized cross-sectional deformations. Gross parameters are para- 
meters which are directly solved for in the displacement formulations 
of ACTION and PLANS. They are not derived quantities like strains, 
stresses, etc. which are derived through spatial differentiation. 

If derived quantities are of importance, one is talking of pointwise 
correlation rather than a global or gross correlation. In such an 
event, the displacement formulations of programs like ACTION or PLANS, 
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even with the inclusion of the effects of cross-sectional felxibility, 
would be inadequate. One will have to then resort to some sort of a 
hybrid or a mixed formulation with the inclusion of large strain ef- 
fects for better response prediction. 
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Figure 1. Thin-Walled Channel Section Column - 

Figure 2. Post-Buckling Response of the Thin-Walled Channel Section 
Column. 
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I. Introduction 


Nonlinear transient analysis of structures has been of increasing 
interest to engineers by virtue of their interest in minimizing human 
and property damage resulting from the catastrophic failure of such 
structures under crash or- seismic conditions. Complexities of the 

' ' I 

structural configuration and its equally complex transient response in 
the presence of material inelasticity make finite element modeling of 
such structures a very natural and plausible recourse. Portions of the 
structures may remain elastic and undergo infinitesimally small deforma- 
tions while other portions may experience finite deformations and mo- 
tions and respond inelas tically under time-varying loads that may lead 
to a complete failure of the structure. If finite strains are to be 
permitted in the model, distinction must be made between undeformed and 
deformed configurations and the concepts of pseudo stresses and con- 
jugate strain measures which have intricate physical interpretations ’ 
must be introduced [1], Furthermore, strictly speaking most elastic- 
plastic theories which hypothesize an additive decomposition of the 
total strain into an elastic and a plastic component lose their validity 
in the large strain domain [2], Because of this, most developers of 
nonlinear analysis codes restrict themselves to a small strain for- 
mulation but permit finite displacements and rotations thereby allowing 
buckling and collapse of the structure to occur. There are some indi- 
cations that this may be adequate for most practical purposes. 

With this hypothesis as its basis, the present discussion focuses 
on the simulation of response of a structure modeled as an assemblage of 
membrane, frame (3-D beam) , stringer elements and rigid links (see 
Figure 1) . The mathematical model is a finite element displacement 
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Figure 1. Finite Element Model of an Aircraft Fuselage Substructure 
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model which consists of discretizing the actual structure by an assem- 
blage of finite elements and approximating the response of each element 
by a finite number of deformation states expressed as linear functions 
of the generalized nodal displacements. 

Two distinct solution approaches exist: (i) the vector approach 

and (ii) the scalar approach. In the former, the mathematical model is 
derived on the basis of the principle of virtual work and reduces to a 
system of nonlinear second-order differential equations in time. In the 
latter approach, a scalar or potential function associated with the 
energy of the model is introduced, minimization of which yields the 
desired equilibrium configuration. In both approaches a temporal finite 
difference scheme is utilized to effectively eliminate time as a var- 
iable. As a result, in the vector approach the equations of motion are 
reduced to a system of nonlinear algebraic equations in the unknown 
nodal parameters of the finite element model ^ In -the scalar 

approach, which is of relevance to this report, the problem is reduced 
to a well known problem in mathematical programming namely the uncon- 
strained minimization of a nonlinear function o'f several variables. 

The scalar approach, has been used for nonlinear analysis by previous' 
investigators [7]-[9]. However, with the possible exception of reference 
[9] .most of the previous work using the energy minimization technique 
has been restricted to static analysis of structures. The algorithm of 
reference [9] had difficulties in'converging to correct solutions because 
of inherent element formulation deficiencies and use of extremely in- 
efficient and expensive finite difference operations for gradients 
besides being restricted to stringer and frame element models only. As a 
result, no meaningful results using the energy minimization approach 
were obtained. The present formulation overcomes such limitations 
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using analytically derived gradients, consistent element formulations 


and the best current variable metric update formula [10] for use in un- 
constrained minimization [11], 
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II. Minimization Technique for Nonlinear Analysis 
a. Formulation Basis 

In this case the problem of response prediction is posed as the 
minimization of a potential function of the unknown nodal parameters of 
the finite element model. For all structural problems with geometric 
and material nonlinearities of the type considered herein such a potential 
function always exists. Although this technique has been hitherto used 
for mainly positive or negative definite systems, other systems which 
fail to be positive or negative definite can be handled by using the 
least squares method or the modified conjugate gradient method with 
preconditioning [12] . In some cases for such systems displacement 
incrementation rather than load incrementation in conjunction with 
conventional unconstrained minimization techniques can also be equally 
effective [13]. 

The minimization scheme as applied to the solution of transient 
nonlinear structural analysis problems consists of minimizing a poten- 
tial function associated with the system for an assumed relationship 
between displacements and time. The displacement- time relation for each 
generalized nodal displacement of a finite element model may be assumed 
of the form [14] 

q e . = 6CAt)\. + (j - SXAt) 2 i| 0i + <At)q 0 . + ^ U-a) 

Y ei ■ Y<4t)ii ei + » - '(') (At)q 0i + q 0i . d-b) 

where q is the i-th generalized nodal displacement at the end of the 
time step and 3 and y are constants. The quantities q and can 
now be expressed in terms of the i-th generalized nodal displacement, 
q^, velocity, q^ and acceleration, at the beginning of the time 
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step and the generalized nodal displacement, q^, at the end of the time 
step. It can be easily verified that the equation of equilibrium for 
an N degree-of-freedom system with lumped masses 


M.q . - F. + ~ - = 0 

i^ex 1 9q . 

ex 


; i = 1,2. . ,N 


correspond to the necessary conditions for the functional 


S = {[ — 

* i=l 20 (At) 


2 q ei 


~ q 0i + >(4t) q 0i + ( 26 


- 1 «0i )q .i ]M i 


- F>. (t+At)q . } + U + C (3) 

x n ex 

to be stationary. In Equation (3), U is the strain energy and C is an 
arbitrary constant. Thus, knowing q^, q^ and'c}^ at time t for any 
given load at time (t+At), the functional S -may be minimized with 
respect to the generalized nodal displacements, q^ (i=l,...N), in order 
to determine the corresponding stable' equilibrium configuration. Thus, 
this scheme satisfies equilibrium at the end of the time step, thereby 
providing an implicit temporal integration scheme. The size of the time 
step is automatically controlled so that the error at half time based 
on interpolated configuration data is less than a prescribed change in 
total energy. In general, the strain energy U will be a nonlinear 
function (at the very least a quadratic) of the generalized nodal dis- 
placements q . Details on the explicit evaluation of U as a function 
of q^ will be touched upon later. 

Of all the available techniques for uncons traxned_ minimization only 
the quasi-Newton or the variable metric methods have been more frequently 
used for structural analysis, because of their higher effectiveness 
[15]. Again, unless one accounts for the sparsity of the Hessian matrix 
and an update scheme which maintains it, one has to almost invariably 
resort to some form of a conjugate gradient technique for problems 
wherein N is an extremely large number. The extension of the minimi- 
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zation techniques to extremely large scale nonlinear structural analysis 
problems is a subject of separate research [16] in itself and is beyond 
the scope of this paper. 

Most algorithms for unconstrained minimization seek a direction of 
travel and the amount of travel in that direction. The manner in which 
these are sought depends upon the sophistication of the particular 
algorithm invoked. Most often the directions of travel are sought in a 
manner which guarantees not only a decrease in the value of the function 
to be minimized at each iteration but also a convergence to the minimum 
in a finite number of iterations (usually N+l for an N dimensional 
space) in the case of quadratic functionals. It is important to note 
that all functionals are very nearly quadratic in the neighborhood of 
the minimum. The present formulation uses the well-known Broyden- 
Fletcher— Goldf arb-Shanno (BFGS) variable metric algorithm which dis- 
penses with the exact line searches while using an Hessian update formula 
which, in the case of a quadratic functional, guarantees a monotonic 
convergence of the eigenvalues of the approximating matrix to the inverse 
Hessian. The iterative scheme is begun with an initial guess which is 
usually the null' vector in the absence of other better estimates. For 
the variable metric or the conjugate gradient methods the required 
gradient of S is evaluated either analytically or by- a finite difference 
operation on S. The use of an analytic gradient results in a substan- 
tial saving in computational effort. This saving is the result of not 
only a cheaper gradient evaluation but most often a faster convergence 
of the solution because of higher accuracy of all computed quantities 
[15]. The i-th component of the gradient of S can be written as 


9S „ . 9U 

t = M.d . - F.+ - — 

9q_. 1 ei x 9q 


ei 


ei 


(4) 
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The term in Eq. (4) requiring significant computational effort is 
3U 

-5 as it embraces the geometric and material nonlinearities. Using 

q ex _ 

half-station central differences this is given by 


9U 

9q . 
ex 


m 


l V?.,. 


k=l 


k VH el’ H e 2 * 


q - + \ Aq ei’ q ei+l * * * W 


ex 


Aq . 
ex 


f U k (q el’ q e 2 
k=l 


*ex 


Aq ei’ q ei+l ‘ * ' 


Aq . 
ex 


(5) 


where Aq^ is a small change in the i-th component and m is the number 
of members or elements which has the i-th degree of freedom in common. 
In evaluating the gradient vector analytically [17], each of its com- 
ponent involves the evaluation of only a single function similar to the 
function for member energy evaluation. Thus, 

3U v t 3W r t ,dW x , 3e 


*, . • l I 

ex k=l v, 


3q 


ex 


dv t- E ^ ^ -^k^q . ^ k dv k 

\ de n ei 


( 6 -a) 


k=l 


where 


W = strain energy density 

2 


. 2 , 2 

(e^ + e - £ £ 


xx 


yy 


,3 2 .1/2 

o + t y ) for a 
xx yy 4 'xy 


e = effective strain 


two dimensional stress state 


„= e for a uniaxial stress state 
xx 


v. = volume of the k-th element 
k 

and for one step incremental loading or unloading 


dW - 
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( 6 -c) 


Equations ( 6 -a) through ' ( 6 -c) imply the use of Hencky’s total strain 
theory along with its assumption that in the strain hardening range the 
inelastic component of the total strain is predominant [18]. This in a 
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way is consistent with the assumption that the total strain can be de- 
composed into an elastic and a plastic part especially in cases where 
the strains are large [2] . According to reference [2-] it is only when 
plastic strains are predominant that such a decomposition is justified. 

The problem at hand could have equally well been formulated using 
the incremental flow theories of plasticity in the strain hardening 
range. The potential function instead of being a function of the total- 
quantities need then be expressed in terms of incremental quanties 
and the minimization technique can still be used [19].’ As a matter 
of fact, it may be conjectured that the performance of the solution 
algorithm will perhaps be significantly improved using such a formula- 
tion even though the material model may then be slightly more complex. 

In any event, it is immediately obvious that a significant reduction in 
computational time will be realized if analytic gradients are used in 
preference to central difference gradients. 

The complexity of the strain energy evaluation for any element is 
determined by its deformation model. This is discussed next, 
b. Deformation model : 

The deformation model of the entire structure is synthesized from 
deformation states of each element of the structure. These states are 
expressed in terms of generalized displacements of the nodes of the 
structure at which the elements interface. The displacement field 
within each element is chosen as a continuously differentiable function 
of the local spatial coordinates and the generalized nodal displace- 
ments. The field maintains interelement continuity of its essential 
derivatives thereby providing a Galerkin model of the system. The local 
generalized nodal displacements of each element are then related to the 
global displacements of the assemblage. These relations, which can be 
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interpreted as transformations of the local coordinate system to the 
global coordinate system, may be linear or nonlinear depending upon 
whether the motions and deformations of the elements are infinitesimal 
or finite. For large rigid body rotations, these transformations are 
accomplished using Euler angles which are linearly independent by virtue 
of the fact that the rotations are performed in a prescribed order. 

There are three kinematic descriptions most commonly used for char- 
acterizing large displacements of finite element models of structures. 
These are: (i) the total Lagrangian formulation wherein the initial 

undeformed configuration is the reference configuration, (ii) the up- 
dated Lagrangian formulation which uses a total Lagrangian formulation 
within each load or time step but updates the reference configuration at 
the end of each step and (iii) the co-rotational or rigid convected 
coordinate formulation which utilizes a coordinate system rigidly at- 
tached to an element and moving with the element. For development of a 
large rotation formulation vital for crashworthiness studies the use of 
the total Lagrangian formulation is unsuitable since most structural 
theories permit only moderately small rotations [20] . The co-rotational 
formulation decomposes the total displacements into a rigid body motion 
component and a strain producing component. Thus, with the restriction 
of small relative rotations within the element, this formulation leads 
to a simplification of the strain-displacement relationship on the 
element level while still permitting arbitrarily large rotations of the 
element. The present deformation model uses the co-rotational or rigid- 
convected formulation for its kinematic description. 

Through appropriate kinematic constraints modeling of massless 
degrees of freedom or of deformation-free rigid links or even the simu- 
lation of contact with an impenetrable, rough plane are easily achieved. 
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Rigid links can be used to simulate either joint eccentrices or rigid 
parts of a structure. In the interest of a truly unconstrained minimi- 
zation Lagrange multipliers or penalty functions are avoided. Rather 
kinematic constraints are formulated as prescribed displacements under 
reactive forces provided by the gradients of the strain energy with 
respect to the corresponding degrees of freedom, 
c. Material model : 

Although closed form analytic expressions for U can be developed 
when the material is elastic the same is not true when elements yield. 
Then the response depends upon the current values of stress components 
and the past history. Von Mises' yield criterion together with Henckey’s 
total strain theory provides a simple means of calculating strain 
energy density distributions throughout an element that has yielded. 
Because total stresses and total strains are no longer linearly related 
recourse must be made to numerical integration (Gaussian or Lobatto) of 
the strain energy density over the volume of the element. 

It is clear then that when an element yields, the complexity of 
the strain energy evaluation increases several times in relation to its 
purely elastic behavior. A number of quadrature points have to be 
assigned over the volume of the element and using the material model 
stresses and strain energy densities have to be evaluated at each of 
these points for known values of strains (Figure 2) . The average strain 
energy density which is simply the weighted sum of these strain energy 
densities then enables the .calculation of the total strain energy. It 
must be noted, however, that the stress-strain history at each of these 
quadrature points, which corresponds to a unique location on an ideal- 
ized effective stress-effective strain curve for the material of the 
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element (Figure 3), must be made available at all times. This places 

4 

highly increased demands on computer storage as inelastic deformations 
progress with time. 

Thus, a frame element which was strictly a uniaxial member in the 
elastic range, typified by its cross-sectional area and moments of 
inertia, requires a full three dimensional characterization in the 
inelastic range. In other words, -frame elements for inelastic analysis 
require a classification based on the different cross-sections. This 
development is restricted to frame elements with thin-walled sections of 
the closed and open (Box, Tube, Elip and E) variety -'a characteristic 
of general aviation aircraft frames. However, in the elastic range the 
development does permit frame elements with arbitrary cross-sections 
characterized by their gross section properties. In the interest of 
simplicity, classical shear flow theory for thin-walled sections is used 
and certain simplifying assumptions regarding torsion, warping and shear 
deformations in the inelastic range are introduced. This is characteristi 
of most nonlinear analyzers mainly because the development of a truly 
three-dimensional frame element for nonlinear inelastic response is a 
formidable task perhaps even more challenging than that of the develop- 
ment of a plate bending or a shell element for the same purpose. In 
fact, in the inelastic range it may be easier to model a thin-walled 
beam of arbitrary cross-section by an assemblage of a large number of 
plate and shell elements thereby permitting a faithful representation of 
very complex effects like restrained warping, torsion, cross-sectional 
distortions, etc. - 

The material is assumed to unload , elastically . For modeling plas- 
ticity under cyclic loading kinematic hardening with an idealized 
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INCREMENTAL DISSIPATIVE 
o- ENERGY DENSITY) 




Bauschinger effect is assumed. Specialized elements like the gap ele- 
ments and stays can be easily modeled simply through an appropriate 
modification of the material model of the conventional elements . 

The details of the total strain energy calculations for the dif- 
ferent element types considered and the transformations relating element 
behavior to global variables is relegated to the appendix. 
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III. Results and Discussion 


' The effectiveness of the minimization technique in solving non- 
linear problems is very much a function of not only the size of the load 
or time step but also the extent and type of the nonlinearity-geometric 
or material and even the type of the temporal discretization scheme used 
which is to say the assumed values of 6 and y in Eq. (I). With this in 
mind, the process of selection of problems for validation was geared 
towards providing an evaluation of the techniques under different types 
of nonlinearities. Problems belonging to three distinct classes namely: 
(i) quasi-static, elastic with geometric nonlinearities, (ii) quasi- 
static, elastic-plastic with geometric nonlinearities and (iii) transient, 
elastic-plastic with geometric nonlinearities were selected. Independent 
solutions or experimental results for these problems were available for 
comparison purposes. 

Figure 4 shows the case of a rod-spring problem wherein the stiff- 
ness of the spring is just enough to prevent a snap- through and provide 
a single-valued load deflection response. Most researchers regard this 
problem as geometrically highly nonlinear. Using stringer elements with 
load steps, as high as 1 lb., the energy minimization solution is indis- 
tinguishable from the easily obtainable exact solution to this problem. 
Higher load steps could have been chosen but caution must be exercised 
with extremely large load steps since the performance (the number of 
minimizations required for convergence) of the minimization algorithm 
may be adversely affected. In other words, the computational effort 
within a load step may increase substantially enough to offset the sav- 
ings accured from fewer load steps. 
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Figure 4. The Rod-Spring Problem 
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LOAD P, kg. 




Figure 5 provides yet another example of such a nonlinearity except 
that in this case the load response curve is no longer single valued but 
is a composite of stable and unstable branches. Using straight-forward 
load incrementation, it is possible to locate only the stable equili- 
brium configurations as indicated in Figure 5. Using displacement 
incrementation, however, the entire load response curve can be easily 
obtained. The response predicted by energy minimization agrees very 
closely with that predicted by the nonlinear analyzer developed by 
Stricklin and Haisler [21] for an identical model of the shallow arch 
using frame elements. 

While both of the previous problems involved only geometric non- 
linearities, Figure 6 presents the case wherein both material and geo- 
metric nonlinearities interact. The experimental prediction of the 
post-buckling, elastic-plastic response of this beam-column with a thin- 
walled channel cross-section was the result of a test carried out by 
Anderson et al [22]. To prevent failure by direct compression, the 
column was tested at an inclination of 5° from the vertical with both 
ends of the column being clamped. A column under these conditions is 
highly imperfection sensitive and hence Anderson et al assume an addi- 
tional 1° offset, as shown in the figure, for the mathematical model 
hoping to simulate the inherent imperfections of the actual column 
tested. Because the response involves a highly unstable branch, dis- 
placement incrementation had to be used in place of load incrementation. 
The response predicted by energy minimization agrees extremely well with 
the experimental prediction and even more so by comparison with that 
predicted using UMVCS-1[23]. 
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Figure 6. Post-Buckling Elastic-Plastic Response of a Thin-Walled Column 
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Figure 7 illustrates the case of the transient' response in the 
presence of both geometric and material nonlinearities. The impulse is 
large enough to cause the entire beam to respond inelastically while 
experiencing moderately large relative rotations. The experimental re- 
sponse for this beam was obtained by Krieg et al [24]. It is immediately 
obvious that the quality of the response prediction is very much a func- 
tion of the values of 6 and y. This is not to say that optimum values 
of 3 and y exist which guarantee optimum fidelity of the response pre- 
diction. As a matter of fact, optimum values of 3 and Y appear to be 
very much problem dependent. Use of 3 - 0.276 and y = 0.55 has been 
recommended by Goudreau and Taylor [25] for smoothing out the high 
frequency oscillations. Again, the response may be significantly 
affected by the use of a consistent mass matrix and with rotatory iner- 
tia and shear deformation effects included. 

With the possible exception of the problem of Figure 8, all the 
previous problems involved only a relatively few degrees of freedom. 
Furthermore, the state of stress in all of these problems was essen- 
tially uniaxial for all practical purposes. Using constant strain 
membrane elements the maximum strain in the vicinity of a notch in the 
direction of loading is determined and compared with the experimental 
results. The agreement between the two predictions is good but. could 
perhaps be improved upon by the use of nonlinear strain displacement 
relationships in the co-rotational coordinate system. 

Next, by way of a reasonably large scale problem consider the drop 
test of a twin engine, low wing Navajo substructure conducted by NASA 
Langley's Impact Dynamics Research Facility under the- auspices of the 
joint NASA-FAA general aviation crash test program. The substructure 
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including the disposition of the occupants and their seats being sym- 

i 

metric about B.L.O.O, a finite element model of only half the structure 
as indicated in Figure 1, suffices. The depth of the substructure floor 
is implicit in the frame elements used for this portion of the substruc- 
ture. Of simulation interest was the occupant chest motion and its 
vertical acceleration at the pelvis location. The occupant was modeled 
by a single lumped mass while the seat was modeled by a set of four 
nonlinear stringer elements whose stress-strain behavior was based on 
previous, independent static crash tests on similar seats. Although, 
the ground plane capability could be exercised for this problem, in the 
interest of simplicity, the aircraft substructure was assumed to be in 
contact with the ground plane at nodes (86) through (9l) while a velo- 
city of 330 inches/sec was imparted to the entire model. Thus one would 
expect correlation in only in the initial phases of the response. 

Figures 9 and 10 provide the correlation between the analysis and 
test results. The chest motion history was obtained from high-speed 
film analysis, performed by NASA Langley. The reduced data for the 
acceleration trace at the pelvis location was obtained by NASA by least- 
square-fit filtering technique. The differences in the two traces on 
the left and right side of the substructure indicate slightly unsyra- 
metrical test conditions. The time of occurrence of the initial peak 
and its magnitude obtained by analysis agrees reasonably well with the 
corresponding test values. For additional model and simulation details 
the interested reader should consult reference, [26]. This reference 
also provides a comparison of the performance of the energy minimization 
technique vis-a-rvis the so-called hybrid technique and another technique 
which utilizes the vector approach. 
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Figure 9. Occupant Chest Motion - Test versus Analysis 
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Figure 10. Occupant Pelvis Vertical Acceleration — Test versus Analysis 



Indeed, one may say that this demonstration of the effectiveness of 
the minimization technique as a tool for nonlinear analysis has been, 
restricted to problems with a relatively few degrees of freedom. For 
such small scale problems the energy minimization technique has been 
shown to be at least comparable to, i‘f not better than, the pseudo force 
technique [16]. Extensions to large scale problems like a full aircraft 
may involve several thousands of degrees of freedom. The state-of-the- 
art in nonlinear transient analysis in general, does not appear to be at 
a point where such large scale problems can be solved efficiently and 
with any high degree of confidence in the simulation fidelity. Likewise 
the effectiveness of the present technique for response prediction of 
such large scale structures remains to be demonstrated. Using precon- 
ditioned conjugate gradient technique or variable metric methods which 
exploit sparsity, it is believed that this is no longer an insurmount- 
able task. 
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VI. Appendix 


Evaluation of the Total Strain Energy 

The scalar approach for the solution of problems of structural 
analysis requires that the strain energy of the system be expressed, 
explicitly or implicitly-, as a function of ‘the global generalized nodal 
displacements of the finite element model. 

From a known vector of the generalized nodal variables in the 
global co-ordinate system, consistent with the prescribed boundary 
conditions, a vector of local generalized variables in the co-rotational 
co-ordinate system of each element is established through transforma- 
tions which are functions of its geometry and its rigid body rotations. 
The assumption of deformation patterns of the element as functions of 
these local generalized nodal variables (interpolating polynomials) 
yields element strains. Recourse- to- the element material model then 
yields the corresponding stresses and strain energy densities at various 
predetermined points (quadrature points) over the extent of the element. 
Barring purely elastic response, a simple weighted summation of these 
quantities over the element volume yields stress resultants and strain 
energies respectively. For purely elastic response these are provided 
by well-known closed form expressions. For the elastic-plastic response 
the strain energy density may be decomposed into an elastic part and an 
incremental-dissipative part thereby providing an estimate of the total 
energy of the system that has been dissipated through inelastic deforma- 
tions. Thus, as shown in Figure 3 for a system with M elements 
M . M M 

U = I U 1 = I U-J + AuJ = J ( / ■ W o 1 dv + / Aw/dv) (A—l) 

i=l i=l e a i-1 v. e v. d 

i 1 - 

where the dissipative energy AU^ is the incremental dissipative energy 
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computed from the previous stress state typied by the point Sq on the 
effective stress-effective strain curve of Figure 3. 

In the following sections expressions for the strain-displacement 
relations are developed for the different elements. 

(i) Stringer Element 

A structural component of uniform cross section which is initially 
straight and which is capable of resisting only axial loads is known as 
a stringer element. 

From Figure (A-l) it can be seen that for assumed nodal displace- 
ments (U , V , W ) and (U , V , W ) of nodes p and q, the change in 

p p p q q q 

length, DL, of the element is given by 

DL = [ (X + U - X - U ) 2 + (Y + V - Y ^V) 2 
q q p p' V q q P P 

+ (Z + W - Z - W ) 2 ] 1/2 - [(X - X ) 2 + (Y - Y ) 2 
qqpp qp qp 

+ (Z - z ) 2 ] 1/2 (A-2) 

which can be simplified to 


DL = L[1 + 


2 (AXAU+AYAV+AzAW) 


2 2 2 
AU +AV +AW 


.1/2 


- 1 « 


(A-3) 


A being the difference operator for q and p end values. Assumption of 
the usual linear interpolation function in the corotational coordinate 
system then yields 

£ - (rj) (A-4) 


(ii) Frame Element (3D Beam) 

A frame element (Figure A-2) is a structural component which is 
initially straight and which undergoes axial, bending and torsional 
deformations resulting from finite displacements and rotations of its 
ends. 
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From Figure (A- 3 ), the displacements of the end q relative to the 


end p can be seen to be 


<5u = (R - R ) - L + (U - U ) 

— — q —p — '■-q ~p 

or in terms of the three components as 

Jsvj = [T] ] 

. 6w / 


(A- 5 ) 


(X - X \ 

q p J 

L 

1 

u - n 
q p J 

Y - Y ) - 

0 

+ [T] 1 

v - v > 

q p( 

1 1 

I 

P 

1 q pi 

\Z - Z ) 

q p 

0 


w - w ) 
q p 


(A- 6) 


where again IL, V_^ and th (i=p or q) denote the global displacements of 

the nodes. The matrix [T]^ can be shown to be [ 9 ] 

[T] = [T_ (cj> , 4. , 4) )][T. (9 , 6 , 6 )] (A- 7 ) 

p 1 T x’ r y v z L 1 xp yp zp 


With 


[T, (ft ,a ,a ) ] = 
l x y z 


c c 

y z 


c s 

y z 


-s 


-c s +s sc cc+sss sc 
xzxyz xzxyz xy 


Ls s +c s c 
x z x y ‘z 


— s c +c s s c c I 
xz- xyz xy 


(A-8) 


c. = cos a. and s. = sin a. for i=x,y and z. Angles <J> , <f> and 4 > are 
1 11 1 & x T y z 


the initial orientation angles described in Figure (A- 2 ) and angles 

9 , 9 and 0 are the rigid body rotations of the end p. In derivin 

xp yp zp 0 J 

Eq. (A- 7 ) Euler angle transformations are implied with the order of the 


rotations being a , a and a . 

& z’ y x 

Similarly, with the restriction of small relative rotations within 

the element, the rotations ijJ , xb and ilj of the end q relative to the 

r x y z n 

end p are 


♦*) 
u = 

V 

With the relative generalized displacements {Su, 6v, 5 w} and 
W x > ipy, ii> z } known the usual deformation patterns of the reference axis 
of the beam element in the co-rotational co-ordinate system are assumed 


IT], 


0 - 0 

xq xp 

0 — 0 

yq ypi 

0 - e . 

zq zp 


(A- 9 ) 




to be 


u(S) = g ~ 

v(5) = f (35 2 - 2 ? 3 )(Sv - z J> ) + ( 5 3 - g 2 )^ 

L S X 2 (A- 10) 

w(S) = f (3g 2 - 2g 3 )(Sw + yji ) - (5 3 - ? 2 H 
i-i s x y 

where g = x/L and y and z are the co-ordinates of the shear center of 

s s 

the cross-section of the beam. The strain of the reference axis can 
then be shown to be 

£=~- nl/f (l~2g) (<5v-z if/ ) + 2(3g-l)i|> 3 

\ (A-ll) 

- ?[f <l-2g)(Sw+y J; ) - 2(35-1)* ] 

i-i o x y 

with n = y/L and £ = z/L. In the above equations it is implicitly as- 
sumed that the lateral displacements and twists are referenced to a 
longitudinal axis through the shear center while the axial displacements 
and rotations are referenced to the centroidal axis. As shorn in ref- 
erence [27] this assumption necessitates the introduction of an addi- 
tional degree of freedom in the axial direction in the interest of equi- 
librium satisfaction in the inelastic range. 

(iii) Membrane Element 

The membrane element of Figure (A-4) is a plane triangular thin 
plate element under constant strain. The element can undergo large 
rigid body motions but its deformation is restricted • to only in-plane 
stretching resulting from finite displacements of its vertices. 

The orientation of the element is uniquely determined by the global 
co-ordinates of its three vertices, p°, q° and r°. The relative dis- 
placements 6u^, 6u^ and dv^ defined in Figure (A-4) can be seen to be 
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(A-12) 


<5u = R ~ R° 

q 

6u^_ = Q cos 3 - Q° cos a 
6v r = Q sin g - Q° sin a 


with 


R° = (X 2 + Y 2 + Z 2 ) lyr2 

rp rp rp 

R = [ (X + U ) 2 + (Y + V ) 2 + (Z +W ) 2 ] 1/2 
rp rp rp rp rp rp 


Q° = <X 2 + Y 2 + Z 2 ) 1/2 

qp qp qp 

Q = r (X + U ) 2 + (Y + V ) 2 + (Z + W ) 2 ] 1/2 
qp qp qp qp qp qp 


(A-13) 


cos a = (X X + Y Y + Z Z )/(Q°R°) 
qp rp qp rp qp rp 


cos 3 = [ (X + U ) (X + U ) + (Y + V ) (Y + V ) + 
qp qp rp rp qp qp rp rp 


(Z + W )(Z + W ) ] / (QR) 

qp qp rp rp 


and typically a^. = - a_. . Next, as in the case of the two previous 

elements, deformation patterns u(x,y) and v(x,y) in the co-rotational 
co-ordinate system when expressed in terms of the local nodal displace- 
ments yield 

<Su <5u Su 

rt v n 

cota]y 

{ 

(A-14) 
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deformation theory as 
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Thus, with the assumption that- the total deformation theory of plas- 
ticity is applicable, the effective strain and effective stress defined 
by Eqs. ( 6— b) and (6-c) yield estimates of the stresses and strain 
energy densities from the material model. It. is obvious that the inte- 
grations over the volume of the element are rendered trivial by virtue 
of the assumption that strains and hence the stresses and strain energy 
densities within the element are constant. 

(iv) Rigid Link 

Rigid link is an element which merely translates and rotates 
without any appreciable deformations. The element is identified by two 
nodes located with reference to the global co-ordinate system. One of 
these two nodes is referred to as the master or primary node, p, with a 
maximum of six independent. degrees of freedom. The motions of the slave 
or secondary node or nodes, q, are determined purely from kinematics by 
setting the left hand side of equation (A-5) to zero. Knowing the 
dependent displacements of the secondary nodes, 

{U} = {U} + [T] T {L} + {R} - {R} (A-16) 

q p p P q 

the contribution, to the total potential energy, of loads applied 
directly at the secondary nodes can thus be determined. 
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T; Introduction 


The prediction of transient linear or nonlinear response of structures is 
almost invariably accomplished by using a temporal finite-difference scheme to 
effectively eliminate time as a variable and reduce the system to a set of 
algebraic equations in the unknown nodal variables of the finite element dis- 
cretization. Finite differencing in time may be either of the explicit or 
implicit type. 

Most nonlinear analyzers linearize response within a time step and use an 
explicit scheme [l]-[3] or an implicit scheme [3]-[4], while a select few do 
not linearize response within the time step and use an implicit scheme [5] or 
an explicit scheme [6]. At the present time there are no clear cut guidelines 
or criteria for the selection of these procedures. Use of hybrid explicit- 
implicit (on-off) schemes have been advocated and is the subject of current and 
future research [7]-[9]. 

For schemes which do not linearize the response within a time or load- 
step, several different techniques for the solution of the nonlinear equations 
may be used. Such techniques have been discussed at great lengths by Bergan' 
[10] and Stricklin et al . [11]. Of particular interest is the technique utili- 
zing the minimization algorithms of mathematical programming. This approach 
has been used successfully for nonlinear structural analyses [5]. [12]-[14]. 

In this case, the problem of finding the solution of the equilibrium equations 
can be equivalently posed as the one corresponding to the minimum value of a 
potential function 

In the past there has been considerable skepticism with regard to the 
effectiveness of the energy minization technique, when compared to other 
methods, such as the pseudo force technique*. It has been also claimed that the 
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extension of the minimization algorithms to large scale problems is virtually 
impossible. The purpose of this paper is to demonstrate by comparison with the 
pseudo force technique that minimization techniques have and can overcome some 
of these objections and that future improvements, in minimization algorithms 
will further improve their efficiency. 

II. Summary of Minimization Techniques for Nonlinear Analysis 

The solution of the equilibrium equations can be equivalently posed as 
the minimization of a potential function. For all structural problems with 
geometric and material nonlinearities of the type considered herein, such a 
potential function always exists. Although this technique has been previously 
used for mainly positive or negative definite systems, other systems can be 
handled by using least squares methods or modified conjugate gradient methods 
with preconditioning. 

The minimization problem as applied to the solution of transient nonlinear 
structural problems (reduction to the static case is obvious) consists of mini- 
mizing a potential function associated with the system for an assumed relation- 
ship between displacements and time. The displacement-time relation for each 
generalized nodal variable of the finite element model may be assumed to be of 
the form [ 16 ]. 

X ei = B(At) 2 X ei + [\ - 3 ) (At) 2 X o1 + (At)* oi + X oi . ( 1 -a) 

' x ei = y(At)x ei + (l - y) (At)x oi + X oi (1-b) 

where X g ^ is the i-th generalized nodal displacement at the end of the time 
step. The 3 and y constants are selected by the analyst to define the inte- 
gration algorithm. It can be easily verified that the equilibrium equation 
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corresponding to the i-tfr degree of freedom of a lumped mass system 


M.X . - F. + 
l ei i 


3U 

3X ei 


= 0 


( 2 ) 


is the stationary condition of the functional 
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i=I 2p(Atr 61 


- (- 


3( At) 


T X . + 
2 oi 


— !_* . + (' 
e(4t) 01 2s 


1)X .)X . ]M. 
oi ei J i 


- F i x ei + U + C (3) 

where C is an arbitrary constant, U is the strain energy, and and F.. 
are the lumped mass and external force, respectively, associated with the i~th 
degree of freedom. 

Thus, once the assumption of the displacement-time relation is made, the 
minimization approach, unlike the incremental stiffness approach, solves the 
actual nonlinear problem within a given load or time step without linearization. 
Consequently iteration at constant load to improve equilibrium or force 
imbalance at the end of a load or time step is not required. 

In the realm of mathematical .programming, the algorithms used for uncon- 
strained minimization can be broadly classified into three distinct classes 
stemming from the level of computational sophistication: (i) the zeroth order 

requiring only function evaluations, (ii) the first order requiring evaluation 
of the gradient as well as the function and (iii) the second order requiring in 
addition a variable metric related to the curvatures of the functional S. Only 
the techniques belonging to the latter two categories have been more frequently 
used for structural analysis, because of their higher effectiveness in compari- 
son with zeroth order techniques [17]. 

For first and second order' methods the' required gradient of S can be 
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evaluated either' by a finite difference operation on S or analyti cal ly . The 
use of an analytic gradient results in a substantial saving in computational 
effort. This saving is the result of not only a cheaper gradient evaluation 
but in most cases, a faster convergence to the solution because of higher 
accuracy of all computed quantities [17]. 

The i-th component of the gradient of S can be written as 


dS —mV - F +■ 

9X . Vei F i + ax . 
ei ei 


(5) 


The term requiring an analytic expression is (|^- ) which can be evaluated 

3X ei 

as 




( 6 ) 


where 


W - the strain energy density for the k-th member 

m - the number of members or elements which have the i-th degree 
in common 

V k - volume of element k 


e = effective strain a 


= ( e 2 +e 2 

\JT xx yy 


. 3 2.1/2 

e xx e yy + 4 Y xy > 


for two dimensional stress state 


e for one dimensional stress state 

XX 


(7-a) 


For one step incremental loading or unloading 


— = a effective stress 1 
de 


, 2 , 2 . o -2x1/2 

- (a +a ~ao + 3x ) 
v xx yy xx yy xy ' 

for two-dimensional stress state 


= a for one dimensional stress state 

XX 


( 7-b) 
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Thus 



(7-c) 


Equations (7-a) through (7-c) imply the use of Henckey's total strain or 
the deformation theory of plasticity. In addition, unloading takes p.lace elas- 
tically with the load path described by the initial elastic portion of the curve 
and for .modeling plasticity under cyclic loading kinematic- hardening is pre- 
scribed. with an idealized Bauschinger effect. 

rv I I 

The term in Eq. (7-c) involves a volume integral which is very 

9X ei 

similar to that required for a member energy evaluation. Hence, each compo- 
nent of the analytic gradient vector involves approximately the same calcula- 
tion effort as one member energy evaluation, two if the node is common to 
two elements, three if the node is common to three elements and so on. Con- 
sequently -a significant reduction in the number of member energy evaluations 
and in CPU time should be realized if analytic gradients are used instead of 
finite difference gradients. This has been vividly demonstrated for problems 
involving different types of non-linearities in reference [17]. 


Ill . Summary of Pseudo Force Technique 
The pseudo force technique, as applied here, uses the initial strain con- 
cept for the treatment of material nonlinear behavior and an incremental up- 
dated Lagrangian formulation for the geometric nonlinear behavior [43], [44]. 
The governing equation at the (n+l)th time step for the transient nonlinear 
analysis of a discretized structural problem is 
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( 8 ) 


= ( 4P „ + 1> + (4 W - + tR n> 

w.ith the following definitions: 
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(9) 


{R^} - .vector of residual forces due to equilibrium imbalance 

{ Au . i > , (Au i } , {ali . i } - increments of displacement, velocity and 

acceleration 


[M] - mass matrix 

{F^> . - vector of internal forces. 


The equations for the static case are obtained by simply neglecting the 
inertia term in Equation (8). For the transient solution a convenient finite 
difference approximation is the central difference integration algorithm. The 
recurrence relations for this temporal operator are 
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where At is the time step. The integration procedure is explicit because 
(Au^ + i } is obtained directly from Eq. (10) using previous information, and it 
is also notable that {Au n+ -|} is obtained from Eq. (8) without factorization 
of the (effective) stiffness matrix in the step-by-step solution. 

Equation (8) is valid for large elastic-plastic deformation provided that 
the appropriate nonlinear terms are included in the strain-displacement rela- 
tions and that the total strain increment can be simply decomposed into elastic 
and plastic components [43]. With an assumption on the plastic strain dis- 
tribution within an element, the effective plastic load increment can be 
expressed as {AQ n+ ^} = [k*] {Ae^-j} where {A& n+ -|} is the increment of nodal 
or element plastic strain and -[k*] is the initial strain stiffness matrix 
used to represent initial strains and to reflect the assumed distribution of 
both total and plastic strains within an element. Incremental plasticity 
relations are used to determine values of stress and plastic strain developed 
throughout the loading history. This technique has the capability of handling 
both strain hardening and ideally plastic behavior. 

IV. Energy Minimization Versus Pseudo Force Techniques 

Structural analysis computer programs differ significantly in their treat- 
ment of elastic-plastic response. Most explicit codes appear to tolerate 
violation of equilibrium in the inelastic range while implicit codes cannot 
get by with such a deficiency. At least two well-known explicit codes, after 
conversion to implicit type using the same element deformation and material 
models, experienced difficulty in converging to a solution and violation of 
equilibrium. For frame elements, the implicit formulation requires an axial 
displacement field in the inelastic range which is an order higher than the 
linear field which is commonly used for a two noded beam-column element in the 



elastic range [19], Furthermore, the number of quadrature points used for in- 
tegration of strain energy densities, stresses, etc. differ from code to code. 
Thus, even if the two finite element models are identical in the undeformed 
state they differ signficantly as -deformation proceeds. This factor has to. be 
taken into account before undertaking an efficiency assessment of the’ two 
procedures. 

An obvious follow-on to the work reported in reference [17] was a 
comparative efficiency evaluation of pseudo force techniques- (incremental 
linearization-explicit) versus energy minimization techniques (nonlinear- 
implicit). For the purposes of Figures 1 through 5, to be discussed below, 
analytic rather than finite difference gradients were employed and Fletcher's 
algorithm [18], more commonly known as the BFGS (Broyden-Fletcher-rGol dfarb- 
Shanno) algorithm, was used for unconstrained minimization of the functional S. 

Figure 1 illustrates the case .of a rod-spring problem involving a single 
degree .of freedom but regarded as a highly nonlinear problem. The technique 
using energy minimization predicted the 'exact' solution using load increments 
as high as 1 lb. (Larger increments could have been tried but were not 
attempted. For quasi-static loading conditions the energy minimization code 
does not provide an automatic selection of load steps guided by error 
tolerances or the like). The pseudo force technique using 1/1 Oth of the load 
increments did not do quite as well. Since it is only a single degree of 
freedom problem a comparison of running times- was not considered meaningful 
although the energy minimization technique was several times faster than the 
pseudo force technique. 
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Next, Figure 2 illustrates the case of the classic snap through of a 
shallow arch. Since the unstable equilibrium branch presents computational 
problems using both techniques, displacement incrementation rather than load 
incrementation was used. The problem using the pseudo force technique w.as run 

1 , ' t ' ^ ^ ’’ , r^, 

on a CDC 6600 machine (taking BO/29-CP0 secs),’ which'.' is known to be about 

■V .. -1 ^ ■ 

2 1/2 times slower than the CYBER 175. The energy minimization method’ took 
9.618 CPU seconds on the CYBER 1 75. A comparison of the equivalent running 
times shows that the pseudo force technique is about five times slower in com- 
parison with the energy minimization technique and the' quality of the response 
prediction does not appear to be as, good. Using strictly load incrementation 
the energy minimization technique using potential energy as the potential 
function can locate only stable equilibrium configurations as indicated in 
Figure 2 by the circles and dotted load path. The pseudo force technique with 
load incrementation becomes singular at the limit load where the tangential 
stiffness vanishes. As a point of reference the square symbols represent a 
first order self-correcting solution from [11]. 

Figure 3 illustrates the post-buckling elastic-plastic response of a thin- 
walled channel cross-section column. Again, displacement rather than load 
incrementation is used because of the unstable equilibrium branch. For the 
range of deformations considered the entire column, built-up from 12 frame 
elements, responded inelastical ly. Figure 3 indicates that the energy minimiza- 
tion technique is comparable to the pseudo force technique .for this problem. 

•Figure 4 provides an example of transient response - in the presence, of 
both geometric and material nonlinearities. The problem consists of a clamped- 
clamped beam of rectangular cross-section subjected to an explosive loading 
over a central region. The experimental data was taken from the literature [411. 
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The impulse is large enough to cause the entire beam to behave inelastically. 

In this case, the amount of nonl inearrty within a time step is very much depen- 
dent not only on the size of the time step but also on the temporal algorithm 
used. That is to say it is very much dependent on the values of 8 and y. 

Thus, it is rather difficult to settle upon a common denominator for comparison 
of the two procedures. In addition, all the previous comments regarding the 
violation of equilibrium in the elastic range and the complications for the 
implicit techniques thereof, hold true for this case as well. With this in 
mind, one would tend to conclude that the efficiencies of the two procedures 
are comparable, although the pseudo force technique appears to have a slight 
edge over the energy minimization technique when the entire structure responds 
inelastically. The particular pseudo force technique code used here for com- 
parison purposes provided the options of using either the central difference 
or Adams predictor^-corrector algorithm. Of these two algorithms the former 
was found to be the more efficient and hence it is used for comparison with the 
energy minimization technique which employed variations of the Newmark-Beta 
algorithm. One should note in passing, however, that the transient response 
appears to be tremendously influenced by the values of 8 and y. This is 
not to say that optimum values of • 8 and y exist which can guarantee optimum 
fidelity of the response prediction. As a matter of fact, optimum values of 
8 and y appear to be very much problem dependent and it may indeed pay-off 
to use a hybrid explicit-implicit (on-off) scheme. 

Finally, Figure 5 provides the case of a structure built-up from two 
dimensional membrane elements. The experimental results for the notched tensile 
specimen are taken from [42]. A comparison of the running times for this 
problem reveals that the time of execution using the pseudo force technique is 
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nearly twice that required by energy minimization. This may be in part due to 
the fact that the pseudo force technique formulates the problem using the 
incremental flow theory of plasticity whereas the energy minimization technique 
uses the simpler deformation theory of plasticity which postulates the exis- 
tence of a strain energy density function of the total stresses and strains as 
in Eq. (7). A comparison of the responses, however, brings out quite vividly 
the effectiveness of both the techniques. It may be remarked in passing that 
the problem at hand could have equally well been formulated using the incre- 
mental flow theory when energy minimization is used as the solution procedure. 
The potential function, in this case, instead of being a function of the total 
quantities need then be'expressed in terms of incremental quantities [20]. 

With the possible exception of the problem of Figure 5 all the other pro- 
blems involve only a relatively few degrees of freedom. However, based on the 
results of Figures 1 through 5 it can be safely concluded that at least for 
small scale problems the energy minimization technique is better suited than 
the pseudo force technique for solving highly nonlinear problems: If anything, 

it is the extent of the inelasticity requiring very costly numerical volume 
integrations which mars its performance. Ignoring the scale of the problem for 
the moment, such a conclusion can have profound implications in the selection 
of a technique for analyzing the highly nonlinear crush response of vehicles • 
wherein inelastic deformations are usually confined to a ‘very small portion of 
the entire structure with a major portion behaving either elastically or like 
a rigid body. 

Results of* Figures 1 through 5 certainly dispel the skepticism with which 
investigators in the past have regarded the effectiveness of the energy minimi- 
zation technique. Needless to say, the solution algorithms of the pseudo 
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force techniques can be considerably improved upon resulting in a much higher 
efficiency. However, by the very repetitive nature of the computations of an 
implicit technique like the ene'rgy minimization technique, savings of as high 
as 50 percent in the computational effort may be- realized simply by an expert 
restructuring of the computer logic alone,. 

It has been also claimed that the extension of the minimization’ algorithms 
to large scale problems is virtually impossible. One would not have challenged 
the veracity of such an assertion back in the sixties, but not any longer. 

Within the past decade mathematicians and computer scientists ha.ve extended the 
scope of the mathematical programming techniques significantly. 

IV. Extension of Minimization Techniques to Large Scale Nonlinear Systems 

Reference [17] concludes that for general "nonlinear structural analysis, 
Fletcher's variable metric method (BFGS) [18] with analytic gradients is the best 
minimization algorithm of those considered therein. This algorithm is initially 
very slow in covergtng to a solution presumably because it uses the null' vector 
as an initial guess for the unknown variables and the identity matrix as an 
approximation to the inverse Hessian matrix [H]-. Thus the solution for the 
first' time or load step is a strong function .of the total number of degrees of 
freedom of the problem. The same is not true, however, of the subsequent time 
or load steps. Having a good approximation to the inverse Hessian [H] and a 
good initial estimate of the variables the number of minimizations required 
for convergence in subsequent steps is only a small fraction of the total num- 
ber of degrees of freedom. Although, Fletcher's variable metric algorithm 
is a very powerful tool its storage requirements (upper or lower half of the 
symmetric matrix [H] requiring n x (n + l-)/2 storage locations for a n 
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degree of freedom system) limits its applicability to small scale problems. 

To alleviate the problems of storage one could fall back upon conjugate gra- 
dient algorithms of Fletcher-Reeves [21] and Powell [22]. But both these 
algorithms have been found to be not nearlyas efficient as the variable metric 
(BFGS) or the Davidon-Fletcher-Powell 's (DFP) algorithms [23] presumably be- 
cause of some ill-conditioning. So an extension of the minimization algorithms 
to large scale problems centers on reducing the storage requirements of the 
second order quasi-Newton methods (BFGS, DFP, etc.) or improving the efficiency 
of the first order conjugate gradient techniques. It is essentially these 
very features which have received the tremendous attention of the mathemati- 
cians and the computer scientists in the past few years. 

The convergence rate of the conventional conjugate gradient method has 
been considerably improved by the use of preconditioning. Briefly, precondi- 
tioning involves the modification of the residuals (components of the gradient 
of the potential function) and working with the pseudo residuals rather than 
the residuals themselves. The modification in designed to reduce the spectral 
condition number of the Hessian of the potential surface and thereby accelerate 
convergence to the minimum. Preconditioned conjugate gradient or generalized 
conjugate gradient methods have been used for both linear and nonlinear problems 
with a great degree of success [24]-[28]. Recently, Axelsson has extended the 
use of generalized conjugate gradient techniques to mixed finite element varia- 
tional formulations [29], Again, using preconditioning Axelsson [29] and 
Widlund [30] have demonstrated the effectiveness of the generalized conjugate 
gradient method in solving problems pertaining to nonpositive definite systems. 
The details on the application of the generalized conjugate gradient technique 
to nonlinear elliptic boundary value problems in irregular regions, as for 
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instance the nonlinear problems of structural analysis using a finite element 
discretization, are found in reference [27]. An additional improvement is an 
attempt to eliminate convergence problems resulting from round-off errors by 
a modification of the linear searches [31]~[32], 

I-n order to be able to use second order quasi -Newton methods for large 
scale problems one has to exploit the fact that the Hessian of the potential 
surface for most fini-te element discretizations is symmetric and banded with 
a very narrow band width. However, the inverse of a banded matrix is not 
necessarily banded and since most variable metric methods in their original 
form required approximations to the inverse of the Hessian the fact that the 
Hessian is banded could not be exploited. Recently however, several investiga- 
tors [33l-[36] have modified the variable metric methods by utilizing anapprox- 
imation to the Hessian rather than its inverse in their march to the minimum 
without destroying the sparsity of the approximation during its updating at 
every step. It has been shown that by doing so the modified method still 
retains the convergence properties of the original method of Broyden [33] which 
did not account for sparsity [37]. This is a major step in the direction of 
extending the minimization techniques for the solution of large scale nonlinear 
problems of structural .analysis. 

Disgressing for a moment-, the differential equations of motion describing 
the nonlinear impulsive response of structures are stiff because of the pre- 
sence of very high and low frequency components .in the initial stages of the 
response. In fact, a boundary layer phenomenon exists such that within the 
boundary layer the use of a standard explicit integration technique will lead 
to stability and round-off problems or in the case of an implicit technique 
will lead to inaccuracies which may destroy, the reliability of the computed 
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transient. Specialized techniques utilizing the concepts from the singular 
perturbation theory have been used to tackle such stiff differential equations 
[37]-[39]. When using minimization techniques to solve such problems the 
boundary layer effect exhibits itself as an ill-conditioned minimization prob- 
lem. Boggs [40] has extended the singular perturbation techniques of Mi ranker 
and others [37]-[39] to modify the conventional variable metric methods so as 
to accelerate the convergence to the minimum value. 

Several other investigators [45]- [47] have advocated the use of self- 
scaling variable metric algorithms. These algorithms use a two-parameter 
family of approximations to the inverse Hessian and determine conditions on 
one of the parameters to improve the condition number of the approximated 
Hessian inverse. The effectiveness of such scaling in conjunction with the 
Hessian updates themselves rather than the inverse updates is a matter that 
needs further investigation. 

Thus, two alternatives are available in extending the minimization tech- 
niques to large scale nonlinear systems, namely: (i) the preconditioned con- 

jugate gradient technique or (ii) the variable metric methods that exploit 
sparsity and utilize singular perturbation theory or scaling to eliminate ill- 
conditioning. The relative performance of these two alternatives is unknown 
at the moment, but it is clear that they both offer the promise of being able 
to solve extremely large scale problems efficiently. Their performance vis a 
vis the pseudo force or other techniques remains to be established. 

VI . Concluding Remarks • 

A comparative efficiency evaluation of the pseudo force technique (incre- 
mental linearization - explicit) versus energy minimization techniques (non- 
linear - implicit) is presented in this paper for the purpose of dispel ing 
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some of the skepticism of the past with regard to the effectiveness of the 
energy minimization technique. The use -of analytic gradients rather than 
finite difference gradients significantly improves the -efficiency of the energy 
minimization technique in the solution* o.f nonlinear structures problems. For 
small scale problems the energy -minimization technique- is better suited than 
the pseudo force technique for solving highly nonlinear problems. 

In the past few years the mathematicians and computer scientists have 
been attacking the problem areas which inhibit the extension of the minimiza- 
tion algorithms to large scale problems. Two alternatives are presently 
available, namely: (i) the preconditioned conjugate gradient technique or 

(ii) the variable metric methods that exploit sparsity .and utilize singular 
perturba'ti on theory or scaling to. eliminate ill-conditioning. 
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The- Rod-Spring Problem 

Snap-Through of a Shallow Arch 

Post-Buckling Elastic-Plastic Response of a 
Thin-Walled Column 

Impulsively Loaded Clamped Beam 

The Notch Specimen 
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